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ABSTRACT 

Digital co-addition of astronomical images is a common technique for increasing signal-to-noisc 
and image depth. A modification of this simple technique has been applied to the detection of 
minor bodies in the Solar System: first stationary objects are removed through the subtraction of 
a high-SN template image, then the sky motion of the Solar System bodies of interest is predicted 
and compensated for by shifting pixels in software prior to the co-addition step. This "shift-and- 
stack" approach has been applied with great success in directed surveys for minor Solar System 
bodies. In these surveys, the shifts have been parameterized in a variety of ways. However, these 
parameterizations have not been optimized and in most cases cannot be effectively applied to data 
sets with long observation arcs due to objects' real trajectories diverging from linear tracks on the 
sky. This paper presents two novel probabilistic approaches for determining a near-optimum set of 
shift-vectors to apply to any image set given a desired region of orbital space to search. The first 
method is designed for short observational arcs, and the second for observational arcs long enough to 
require non-linear shift- vectors. Using these techniques and other optimizations, we derive optimized 
grids for previous surveys that have used "shift-and-stack" approaches to illustrate the improvements 
that can be made with our method, and at the same time derive new limits on the range of orbital 
parameters these surveys searched. We conclude with a simulation of a future applications for this 
approach with LSST, and show that combining multiple nights of data from such next-generation 
facilities is within the realm of computational feasibility. 

Subject headings: Solar System; Data Analysis and Techniques 



1. INTRODUCTION 

Detecting faint objects in the Solar System is a techni- 
cally challenging problem. The approaches taken in the 
past can be broken into two basic categories: inter-frame 
source detection and linking, and the "shift-and-stack" 
approach. Inter-frame detection, where transient sources 
are detected in individual exposures and motion is linked 
between multiple exposures, is widely employed and fa- 
vored in planned proposed surveys like Pan-STARRS 
(Denneau et al. 2007) and LSST (Axelrod et al. 2009) 
due to its simplicity and robustness. This method, how- 
ever, does not exploit imaging data to the fullest extent 
possible, since detections are subject to the completeness 
limit of individual frames. 

Since Solar System objects are moving across the sky, 
their images trail and the noise contribution from the 
sky background increases. Trailing losses create an up- 
per limit to individual frame exposure times for efficient 
detection of Solar System objects. It is advantageous 
to combine multiple exposures, each short enough that 
trailing effects are negligible. Co-addition is routinely 
done to improve the SN of stationary sources and to re- 
move cosmic ray events, but in order to apply it to the 
detection of moving objects several other steps must be 
performed. First, contamination from stationary sources 
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can be removed by subtracting a high-SN template im- 
age from each image in the set. Next, the sky motion of 
the objects of interest must be predicted and this motion 
compensated for by shifting each image's pixels in soft- 
ware. When the images are then combined only the flux 
from sources that moved at the predicted rates will be 
constructively added. 

This method first appeared in literature in Tyson et 
al. (1992), though details of the application were lim- 
ited. An early detailed description of this technique was 
presented by Cochran et al. (1995). A set of images from 
the Hubble Space Telescope (HST) Wide Field Planetary 
Camera 2 (WFPC2) were combined to allow statistical 
detection of very faint TNOs. The sky motion of the 
sources was parameterized as angular rates (9) and an- 
gles on the sky (4>). Gladman et al. (1997) used the same 
technique in order to combine a series of images from the 
Palomar 5m telescope in the hope of constraining the size 
distribution of Kuiper Belt Objects (KBOs) by discov- 
ering objects with smaller radii than had been detected 
in most previous wide-area surveys that employed inter- 
frame detection. This survey was followed by further 
searches from Keck (Luu & Jewitt 1998; Chaing & Brown 
1999), CTIO (Allen et al. 2001 and 2002; Fraser et al. 
2008), Palomar (Gladman et al. 1998), VLT and CFHT 
(Gladman et al. 2001), which employed the same tech- 
nique to varying degrees of success. Recently, Fraser & 
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Kavelaars (2009) and Fuentes, George, & Holman (2009) 
both used similar techniques to search Subaru data for 
KBOs as faint as R ~ 27. 

All of these surveys had observational arcs less than 
two days in length. Over periods this short, the mo- 
tion of outer Solar System sources can be well approxi- 
mated (with respect to ground-based seeing) by the lin- 
ear trajectories that the method these surveys imple- 
mented implicitly assumes. For longer arcs or higher- 
resolution data, however, the approximation of a linear 
trajectory is no longer adequate. Bernstein et al. (2004) 
employed a more advanced approach, which they called 
"digital tracking," to extremely high-resolution HST Ad- 
vanced Camera for Surveys (ACS) data. Even in arcs on 
the order of a day in length, nonlinear components of 
motion were non-negligible relative to the ACS PSF. In 
order to stack frames such that distant moving sources 
add constructively in this regime it becomes necessary 
to generate nonlinear shift-vectors via full ephemerides 
from orbits of interest. They stacked arcs up to 24 
hours in length using shift-vectors parameterized as a 
grid over three parameters: two linear components of 
motion and heliocentric distance d. By sampling densely 
over this grid, they were able to recover orbits with 25 
AU < d < oo and i < 45°, and detect objects as faint as 
R ~ 28.5. However, this dense sampling, somewhat un- 
constrained discovery volume, and extremely small PSF 
required ~ 7 x 10 5 shift- vectors, which translated into 
~ 10 pixels to be searched. Bernstein et al. found 
that stacking more than 24 hours of their data became 
computationally prohibitive. 

In the age of ground-based gigapixel CCD mosaic im- 
agers, how is data obtained over long arcs going to be 
efficiently exploited to search for faint sources? In the 
interest of exploring massive datasets with long observa- 
tion arcs to as deep a limit as possible, we have created a 
simple probabilistic method for generating an optimized 
set of shift-vectors for any imaging survey and any in- 
teresting volume of orbital element space. By generating 
shift- vectors in a probabilistic way, exploring an explic- 
itly demarcated region of orbital parameter space be- 
comes straightforward. The targeted region of orbital pa- 
rameter space is searched with uniform sensitivity, sim- 
plifying the accurate characterization of any detected 
sample of objects. This method also avoids searching 
non-physical or uninteresting parameter space, thereby 
maximizing scientific return for minimum computational 
cost. 

In Section |2.2[ we define the maximum tracking er- 
ror, the basic parameter that defines the density of a 
grid. Section [2] describes previous analytical methods for 
search-grid parameterization, our probabilistic approach, 
and a grid-tree structure for optimizing multiple-night 
arcs. Section [3] compares the efficiency of our method to 
those used in previous surveys, and in Section|4]we show 
an example application of this method to the kind of 
data that may be acquired by the Large Synoptic Survey 
Telescope (LSST). 

2. GENERATING SEARCH GRIDS 

2.1. Estimates of On- Sky Motion 

Previous on-ecliptic surveys (eg., Fraser & Kavelaars 
2009) have used simple analytical estimates for the ap- 



parent on-sky rates of motion for distant Solar System 
objects in order to generate shift rates. For an on-ecliptic 
observation of a field /3 degrees from opposition, the 
range of on-sky angular rates of motion 9 and angle from 
the ecliptic <f) for a distant object at heliocentric distance 
d and geocentric distance A with inclination i and ec- 
centricity e can be approximated by the vector addition 
of the reflex motion of the object due to Earth's motion 
(inversely proportional to A) and the intrinsic on-sky or- 
bital motion of the object by Kepler's law (proportional 
tod- 3 / 2 ): 
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where Vd — ^/lzte, representing the fraction of the 
mean angular velocity of an orbit at pericenter (+e) or 
apocenter (— e), compared to a circular orbit at heliocen- 
tric distance d. 

This motion can also be parameterized as two compo- 
nents of angular rates of motion, one parallel and one 
perpendicular to the ecliptic. Re-arranging Eqn. [l] into 
parallel and perpendicular components, we find: 

0„ ~ 148(cos(/3)(A)- 1 - v d d~i cos(i)) "hr- 1 (3) 

9\ ~ U8v d d~l sin(i) "hr- 1 (4) 

The 9 and <fi parameterization is c onvenient for visual- 
ization purposes, but in Section 
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we show that the 9 U 
and 9± parameterization is more appropriate for gener- 
ating a search grid. 

The angular rate of motion is lowest (for a given d and 
A ~ d) for i = 0° orbits at pericenter, since v d is at 
its highest and the parallax and orbital motion vectors 
are anti-alignedj]] The highest angular rate of motion 
at opposition is reached in two cases, depending on the 
maximum inclination and eccentricity considered (i m ax 
and e max ): 
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Given this approximation, the highest <f> any orbit can 
achieve is a strong function of d. If we consider an orbit 

1 This assumes that the object's distance and eccentricity are 
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not such that that v^jd^ > - in other words, the object is not 
overtaking Earth at opposition. If this were the case, 9 is lowest 
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with e — > 1 at pericenter with i = 90°, and approximate 
A ~ d, then by Eqn. [2] it can be shown that: 
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2.2. Grid Geometry and Maximum Tracking Error e max 

In using digital tracking, the motions of real sources 
are approximated by a grid or some distribution of shift- 
or motion- vectors. In order to quantize the effectiveness 
of this strategy, we define a maximum tracking error e to 
be the largest possible error tolerated between any real 
object's motion and the nearest predicted motion vector. 
To estimate its effect, we compare the final signal-to- 
noise ratio S/N' for a faint circular source with flux / 
and area Aq before and after a linear tracking error e 
is added in an image with exposure time t. The linear 
tracking error e adds an additional rectangular region 
to the area of the source, with width 2R (where R is the 
aperture radius applied to the initial circular source) and 
length e. This rectangle has area 2Re, and the total area 
of the blurred source is now: 



A' = ttR 2 + 2Re = A n [l + —e 
The final signal-to-noise S/N' is given by: 
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So, for a search-grid with maximum tracking error e max , 
the worst possible signal-to-noise degradation is a factor 
of F = (I+^e)~2. Given a gaussian point-spread 
function with a full- width at half-maximum of T, the 
radius at which the signal-to-noise of a sky-dominated 
source is maximized is R ~ 0.68I\ If we define this as 
our initial aperture radius, then the maximum allowable 
tracking error given a desired limit on F becomes 



e ma x - -0.68T(F- Z - 1). 



(7) 



Previous surveys have used widely varying values for 
f-max with respect to the typical seeing of that survey, 
ranging from ~ 0.1T — 3r. Surveys that searched data 
by eye were limited in the number of shift- vectors they 
could apply, and leveraged the eye's robust noise discrim- 
ination to compensate for the loss in signal-to-noise due 
to using large e max (eg., Fraser & Kavelaars 2009). Sur- 
veys using automated detection pipelines have decreased 
£max to limit their signal-to-noise loss, and leveraged 
massive computational resources to compensate for the 
increased number of required shift- vectors (eg., Fuentes 
et al. 2009). 

In order to evaluate the maximum tracking error of any 
search-grid of shift- vectors applied to a given data set, we 
determine the distribution of shifts from the initial image 
to the final image, where the shift for orbit k in image i 
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The maximum separation any real orbit's motion from 
one of these offsets determines the maximum tracking 
error of the search-grid. After time t has elapsed, the 
on-sky angular separation SI of two objects starting from 
the same initial position with identical angular rates of 
motion 9 and on-sky angles of motion separated by small 
angle d<j) is given by 



SI = 26t sin 



(9) 



The maximum tracking error between a search-grid of 
shift- vectors and a real distribution of orbits depends on 
the geometry of the search-grid. Consider the grid spac- 
ing of a geometrically regular search-grid in the plane 
defined by the total changes in RA and DEC made by 
moving sources during the observations, which we will 
refer to as the (da, dS) plane. If the points of the search- 
grid are arrayed at the vertices of identical adjoining geo- 
metric cells on the (da, dS) plane, the point most distant 
from any point on the search-grid is the point at the cen- 
ter of each cell. The distance from any vertex of a cell 
to the center of that cell defines the maximum tracking 
error of that search-grid. 

For points arrayed on the vertices of a grid of adjoining 
rectangles with sidelengths dA and dB, the maximum 
tracking error is e max = 'dA 2 + dB 2 . In the case there 
dA = dB, this becomes e max = \[2dA/2. This geometry 
well approximates the search-grids used by most previous 
surveys. The left panel of Figure [I] illustrates a grid with 
this geometry. 

If a fixed (dd,d(f>) grid spacing is adopted such that 
a maximum tracking error e max is satisfied after time 
t for objects with angular rate of motion 9 = 9 , ob- 
jects with 9 > 9q will see increasing tracking errors. 
To illustrate, consider a roughly rectangular grid where 
dA ~ d9t and dB ~ 29tsm(d(f>/2). Inserting these val- 
ues into our definition of e max for a rectangular grid, we 
see that e max oc 9. A fixed (d9,dcf>) grid parameterization 
over-samples motions with low 9, and under-samples mo- 
tions with high 9. 

The over- and under-sampling problem can be reme- 
died by parameterizing the sky motion as two compo- 
nents of angular rates of motion, one parallel and one 
perpendicular to the ecliptic (9 n and 6±). By adopting 
a fixed spacing in (6 n , 9±) along with angles from the 
ecliptic (4> m i n , 4>max) between which motions are con- 
fined, a uniform e max can be maintained for all angular 
rates of motion. This is the approach that was adopted 
by Fuentes et al. (2009). 

The packing efficiency of a search-grid can be improved 
over the regular rectangular case. If the area of (da, dS) 
to be searched by a grid is much larger than the unit 
of area searched by an individual grid point such that 
we can largely ignore boundary conditions, the proof by 
Gauss (1831) regarding the most efficient regular pack- 
ing of circles on a plane holds. In this case, defining 



4 



Parker & Kavelaars 



the points of the search-grid as the vertices of a grid of 
adjoining equilateral triangles will produce the highest 
packing efficiency. Considering a grid of equilateral tri- 
angles with sidelength dA, the maximum tracking error 

is Cmax = dA/y/3 (where the effective dB spacing be- 
comes 3e max /2). The right panel of Figure [I] illustrates 
a grid with this geometry. Both our short- and long-arc 
solutions described in the next section produce a grid 
with approximately this geometry, and typically require 
> 20% fewer grid points than a similarly-defined rectan- 
gular grid with the same e max . 

Regardless of the grid geometry, each cell is basically 
covering a small area of the final (da, dS) plane. The 
final area of the (da, dS) plane covered by real objects 
is approximately oc t 2 (Equations |3 and [4[ neglecting ac- 
celeration), while the area covereaby a single grid point 
is oc e^j aa . The total number of grid points N required 
to search a given population scales as t 2 /e 



2 

max ' 



2.3. Probabilistic Methods for Generating Shift-Vectors 

The methods for generating shift-vectors described 
above are cumbersome to analytically generalize over any 
region of the sky, arbitrary regions of orbital parame- 
ter space, and arbitrarily long observational baselines as 
they still suffer from the problem of non-linear compo- 
nents of motion that become significant as t grows. In 
this regime it becomes necessary to add additional di- 
mensions to the search grid (9 and <fi, for example), and 
selecting an optimized set of grid-points becomes a chal- 
lenge. Instead of pursuing an analytical solution, in the 
following section we outline a simple probabilistic ap- 
proach for defining an optimum set of shift-vectors for 
any region of the sky and over any length of observa- 
tional baseline at any solar elongation. 

Instead of approximating the minimum and maximum 
angular rates of motion and angle on the sky for an ad- 
hoc subset of the orbits of interest, our method takes the 
following approach: First, generate a very large sample of 
synthetic orbits that fill (in a characterized way) the or- 
bital parameter space of interest. Second, use ephemeris 
software to determine the (a, S) shift between frames for 
every synthetic orbit in every imaged epoch. Finally, 
search for the optimum (smallest) subset of orbits whose 
motions accurately represent the motions of the entire 
set of synthetic orbits within a maximum tracking error 
tolerance. The set of (a, 6) shifts (translated into (x, y) 
pixel shifts) from this representative set of orbits is our 
optimum set of shift- vectors. A more detailed description 
of the method follows. 

2.3.1. Initial sample generation 

To create our initial sample, we use a randomly gen- 
erated set of synthetic orbits that fill a region of or- 
bital parameter space of interest such that the orbital 
ephemerides intersect the imaged field during the times 
of observation. The approach to filling orbital parameter 
space is subject to some consideration, as any orbital- 
space biases may translate into a skewed on-sky motion 
distribution, resulting in a biased selection of "optimum" 
orbits. Over small ranges of d, a uniform sampling is ad- 
equate, but for larger ranges we sample uniformly with 
respect to d~ c , where the index c = 2 is appropriate to 
generate a uniform distribution in 9 for observations near 



opposition. Once we have selected a given d, we select 
a, and e from a uniform distribution between the mini- 
mum and maximum values set for those parameters, in 
the following order: first select e, then a such that it falls 
within its limits and also a > d(l — e). 

In cases where a population has a particular orbital pa- 
rameter that is well-constrained, it is preferable to gen- 
erate this parameter first (eg., a for a mean-motion res- 
onance), and then force the other parameters onto as 
uniform a grid as possible. In the case of populations 
in mean-motion resonances, a is pinned at the resonant 
value a mmr , then d is selected from a uniform distri- 
bution between a mmr (l - e max ) and a mmr (l + e max )- 
Finally, e is selected such that a mmr (l — e) < d < 
immr(l+e). In either case, a, e, and d define two possible 
values for M , which are given equal probability. 

Sampling inclination from the distribution of an 
isotropic sphere (probability of i oc sin(i)) may seem ap- 
propriate so as not to bias the on-sky motion distribution 
to that of low inclination orbits, but because the com- 
ponent of motion perpendicular to the ecliptic is roughly 
oc sin(z), and because the small area of any imaged field 
with respect to the entire sky places stronger limits on 
the phase-space volume available to high-inclination ob- 
jects than that available to low-inclination objects, we 
contend that is preferable to sample i from a uniform 
distribution. The minimum inclination orbit possible to 
be observed is i min = arccos(? mi „)(l — g), where l min is 
the minimum ecliptic latitude of the field. In the case 
where retrograde orbits are of interest, duplicate the set 
of inclinations with i r = 180 — i. 

The remaining orbital elements to be generated are 
the node and the argument of perihelion. These can be 
rotated to force the object to fall on the imaged field dur- 
ing the dates of observation, and should fill the parame- 
ter space set by these constraints smoothly. M and the 
argument of perihelion are coupled for objects in mean 
motion resonances, but for the purposes presented here 
it is convenient to simplify the treatment by treating the 
two as independent. 

Depending on the size of the orbital region of interest 
and the length of the observational arc, the number of 
synthetic orbits required to produce a statistically ade- 
quate sample varies. The basic requirement is that the 
density of points in the final (da, dS) plane is such that 
the max spacing between any two nearest neighbors in 
the initial sample is -C e max , ensuring that there are no 
empty bins in (da, dd) given circular bins with radius 
(■max- In trials, we found that this required the initial 
sample size to range from several thousand to hundreds 
of thousands of synthetic orbits. 

Once a set of synthetic orbits has been created, the 
shift-vectors for each orbit can be determined. We de- 
fine the shift-vector for orbit k in image i to be a point 
in the (dai, ddi) as given by Vfc^ in Eqn. ^1 The list of 
shift-vectors [vfc,o...n] is stored as a list linked to the orbit 
k for which they were generated. We then compare each 
list of shift- vectors to determine the optimum set, which 
is the smallest subset which have motions that accurately 
represent the whole set within some spatial tolerance in 
the image plane. We set this tolerance by specifying that 
all synthetic orbits in the initial sample must be matched 
to the motions of at least one orbit in the optimum sam- 
ple to within our maximum tracking error tolerance e max 
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Fig. 1. — Comparison of grid geometries. Left panel: Packing of 
hexagonal lattice, with scales illustrated. 

in every imaged epoch. In other words, every synthetic 
orbit k must be matched to at least one optimum orbit 
h in every image i by |v M - v M | < e max . 

2.3.2. Short- arc solution 

For observational baselines over which any nonlinear 
component of motion is negligible compared to e, the 
largest offset between any orbit and the nearest predicted 
motion vector will always occur in the final image n. To 
generate shift-vectors in this regime, we use the following 
approach: 

1. Propagate all synthetic orbits in the initial sample 
forward to the final image (n) plane and record 
the final shift-vector v^. „ for each orbit k. Find 
the median of the distribution of shift-vectors and 
set the corresponding point in (da n , d5 n ) to be the 
center of a geometrically-optimum triangular grid 
with sidelength dA = V3e max that extends over 
a region of (da n , dS n ) significantly larger than the 
that occupied by any synthetic orbits in the initial 
sample. 

2. From this initial large grid, select only those grid 
points that have one or more synthetic orbits whose 
final shift-vector lies within e of it in (da n , d5 n ). 
Record the list of orbits that are matched to each 
grid point. 

3. After the initial grid has been cut down to just 
those grid points that are matched to at least one 
synthetic orbit, find the smallest subset that are 
matched to unique synthetic orbits: First, select 
the grid point that is matched to the largest num- 
ber of synthetic orbits, set it and its list of orbits 
aside, then remove the references to each synthetic 
orbit recovered by that grid point from all other 
grid points. Repeat this process for the remaining 
grid points, until all the original synthetic orbits 
are accounted for in the set-aside list. All set-aside 
grid points are now identified as the optimum set 
for recovering the initial sample of synthetic orbits, 
given this grid geometry and orientation. 

4. Iterate Steps 1-3 several times, each time rotating 
the initial grid around its central point by some 
small angle d^i (up to a total rotation of 60°). 
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a square lattice, with scales illustrated. Right panel: packing of a 

Record the number of optimum grid points re- 
quired for each orientation, and select the orien- 
tation that requires the fewest. As a small refine- 
ment, if any of the optimum grid points arc found 
to lie with centers outside the distribution of syn- 
thetic orbits in the (da n , dS n ) plane, we shift its 
center to the nearest (da n , dd n ) point from a syn- 
thetic orbit. After this refinement, the set of opti- 
mum grid points from this orientation arc defined 
as our final optimum set. 

These grid points are points in the final (da, dS) plane. 
We can back out the implied set of optimum (9 a , 6g) 
rates by simply dividing the respective components of 
each grid point by the total observational baseline tb- 

The rotation in step 4 is included to address boundary 
conditions. If the region of (da, dd) searched by a grid is 
not ^> e max , then there will be a preferred orientation of 
the grid. However, if e max is small relative to the total 
region, this step is unnecessary. 

If, instead of trimming the search-grid to as small a 
number of points as possible, it is preferable to create 
a "safe" grid that over-searches the boundaries in rate- 
space by some comfortable margin, then modify step 2 
to select grid points that lie within a larger distance N x 
£ max of any final shift-vector. Step 3 must be skipped in 
this case. This maintains the grid spacing and geometry, 
but adds a "buffer" region around each synthetic orbit. 

2.3.3. Long-arc solution 

In the case that the observational baseline tb is long 
enough such that non-linear components of motion are 
non-negligible with respect to e, a different treatment 
is required. Instead of presuming that the largest offset 
between the motions of a given orbit and the nearest grid 
point will occur in the final frame n, we must instead 
ensure that for every synthetic orbit there is one non- 
linear shift-vector which matches its motions within e in 
every frame. While nonlinear components of motion are 
non-negligible, the dominant factor driving the number 
of shift-vectors required is still the linear components of 
motion. 

1. After creating an initial sample of synthetic orbits, 
generate a search grid via the short-arc (linear) so- 
lution for the observations. When generating the 
search grid, set the effective maximum tracking er- 
ror to be slightly smaller than the desired final 
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maximum tracking error. This is to take into ac- 
count that the final grid will be selected from a 
discrete distribution; namely, the motions of the 
synthetic orbits in the initial sample. 

2. For each grid point p on the optimum linear search 
grid, take the list l p of orbits uniquely matched to 
that grid point and consider their shift-vectors in 
every imaged epoch. For each orbit k in list l p , find 
all the other orbits within l p that have motions that 
are matched to the motion of k within the desired 
maximum tracking error e max in every frame. 

3. Select the single synthetic orbit fci with motions 
matched to the largest number of unique synthetic 
orbits in l p and set it aside as an optimum orbit. 

4. If some of the orbits linked to the grid point p are 
not recovered by the first optimum orbit k\ , remove 
all orbits recovered by the orbit k\ from considera- 
tion, then repeat Step 3 to select a second optimum 
orbit fc 2 . 

5. Repeat Step 4 until all the orbits linked to the grid 
point p are accounted for by one or more optimum 
orbits. 

6. Repeat Steps 2 — 5 for each grid point in the initial 
solution. 



£j = Cmax/2. Thus, while this 2- level case requires fewer 
pixel additions, it results in significantly more pixels to 
search than in the 1-level case. 

While we have not implemented this method, we can 
estimate its results. Since the number of grid points N is 
approximately proportional to t 2 /e 2 , and the 2-level case 
requires e, = e max /2, we can estimate that this method 
will require four times as many shift-vectors per unit time 
interval as a 1-levcl grid. If we consider n nights of data, 
where the length of a night is equal to the length of a day 
such that the total observational baseline is tb — (2n — 1) 
in units of the length of a single night. Thus, if the 1-level 
grid requires Ni shift-vectors for the full observational 
timeline, we can estimate that the 2-level grid for the 
full length of time will require N 2 j u U — 4Aq. The 2-level 
method will require ~ 4 times more pixels be searched 
in the final stacks than the 1-layer method. However, to 
determine the improvement in pixel additions, we need 
to determine the size of each nightly grid: 



2,nightly 



2, full 



4Aq 



(2n- l) 2 



(10) 



For surveys where the number of observations taken 
per night N obs > N 2 j u u/N 2 . nighH y = (2n - l) 2 such 
that the number of pixel additions is dominated by the 
nightly stacks, we can estimate the total number of pixel 
additions required by: 



The shift-vectors of the final optimum set of orbits are 
then identified as the optimum set of non-linear shift- 
vectors to apply to the given data set. 

2.3.4. Grid Tree: Optimization for multi-night arcs 

Other significant optimizations can be made to reduce 
the total number of pixel additions required to search 
a given data set. One particular optimization that our 
probabilistic method lends itself to is the use of a multi- 
level grid tree described by Allen (2002). This approach 
recognizes that for multi-night stacks, there are essen- 
tially two distinct timescales: the unit time periods over 
which images are obtained (length of a single night), and 
the total observational baseline. Allen (2002) points out 
that a significant reduction in total number of required 
shift-vectors can be obtained if, instead of combining all 
images obtained over several nights with one search-grid, 
images are combined on a per-night basis (with a single- 
night grid), then these stacks are combined with a sec- 
ond search grid to remove the motion that would have 
occurred over the entire period. If all else remains equal, 
this tree-of-grids structure results in fewer overall pixel 
additions, as redundant combinations of images taken on 
a single night are only performed once. Since this method 
relies on combining images that are the combination of 
other images, it can only be applied for combination al- 
gorithms that can nest: for example, weighted averaging 
and co-adding are usable, but not medians. 

It is important to note that in combining a tree of grids, 
the resulting maximum tracking error e max goes as the 
sum of each level's tracking errors X) e i- As sucn > each 
level's tracking error must be limited to ej = e max /N, 
where N is the total number of levels. In the case 
described here, the tree has two levels, and so the ef- 
fective tracking error in each level must be limited to 



^2, add — N2,nightlyNpi xe l s N b s n 

_ ANjNpteetsNobsn (11) 

~ (2n-l) 2 

Whereas for the 1-level grid, the number of pixel addi- 
tions required is N 1Md = N\N pixels N oba n. The resulting 
improvement in the number of required pixel additions 
is roughly 

NiMd (2n-l) 2 , , 

N 2 ,add 4 1 > 

So for a two-night arc, we estimate a factor of ~ 2.25 
improvement in the number of pixel additions, while for 
a three-night arc this increases to ~ 6.25. 

3. COMPARISON OF PROBABILISTIC SOLUTION TO 
PREVIOUS SURVEY GRIDS 

In the following section, we will apply our probabilistic 
shift-vector generation method to previous surveys' ob- 
servations, and compare the results to the methods used 
by the authors of the surveys. In order to determine the 
limits of parameter space searched by each survey, we 
generate our initial sample of synthetic orbits over very 
large ranges of heliocentric distance (20 — 500 AU), incli- 
nation (0° - 180°), and eccentricity (0 - 0.999). The sky 
motion of each synthetic orbit generated is then tested 
to ensure that it is within the (0,<j>) or (#,,, 8±) ranges 
searched by the authors; as such, we only perform this 
characterization for surveys for which we can accurately 
reproduce the original search-grid from the literature. 

Three surveys are selected for comparison: Fraser et al. 
2008, Fraser & Kavelaars 2009, and Fuentes et al. 2009. 
These surveys are selected because they cover relatively 
large areas (0.25 — 3 square degrees) to relatively faint 
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Table 1: Survey Parameters 



Survey & Method 


Baseline & Tracking Error 


Grid Limits 


Spacing 


^ qrid 


Fraser et al. 2008 


% ~ 4 hours 


V'Ahr- 1 < 


< 4".lhr- 1 


d9 = 0".7hr- 1 


25 


(searched by eye) 


e = 1".6, F ~ 0.57 a 


-10° < 


< +10° 


d(j> = 5° 




Fraser & Kavelaars 2009 


t\y ~ 4 hours 


0".4?tr- 1 < 


9 < ^'.bhr- 1 


dd = 0".21hr- 1 


95 


(searched by eye) 


e = 1".25, F ~ 0.61 a 


-15° < 


4> < +15° 


d(f> = 7.5° 




Fuentes et al. 2009 


ti, ~ 8.5 hours 


0".7hr~ 1 < I 


% < 5'M/ir" 1 


d6„ = 0".lhr~ 1 


736 


(automated detection) 


e = 0".6, F ~ 0.76" 


-l".4hr- 1 < 6 


ij_ < +l".4hr- 1 


d9 ± = O'M/ir" 1 








-15° < 


4> < +15° 







A.<J \ </> \ T-LO 

Maximum S/TV degradation factor due to tracking error, from Eqn. |7| S/TV' = F X S/N 



limits (i? ~ 26 — 27), and because they contain very com- 
plete information regarding their respective search grids 
and targeted orbital parameter space, which are outlined 
in Table 1. The comparisons between the original and 
optimized grids and the heliocentric distance range our 
analysis indicates each survey was sensitive to is listed in 
Table 2. 

3.1. Fraser et al. 2008 

Fraser et al. 2008 searched ~ 3 square degrees (taken 
over baselines ranging from 4 to 8 hours) for TNOs using 
fixed grids of rates and angles, searching the final stacks 
by eye. The images were acquired from several facilities, 
and detailed grid information is only supplied for the 
MEGAPrime observations. These observations spanned 
~ 4 hours, and the search grid applied (described in Ta- 
ble 1) required 25 rates and angles be searched. The 
adopted search-grid spacing was fixed in d6 and d(j>, re- 
sulting in a maximum tracking error as a function of 9 as 
described in Section |2.1| We estimate that the resulting 
maximum tracking error was ~ 1".6 ~ 2.2r, resulting 
in a maximum S/N degradation factor of F ~ 0.57. As 
the authors reported no sensitivity loss as a function of 
rate of motion, we adopt this as our uniform maximum 
tracking error. The authors state that their selection of 
rates and angles were designed to detect objects on cir- 
cular orbits with heliocentric distances from ~ 25 — 100 
AU with inclinations as high as 70°. 

Due to the short arc and large tracking error of this 
survey, our optimum grid (with 19 shift-vectors) is not 
radically more efficient (~25%) than the original grid of 
25 shift- vectors. The original and optimum grids, along 
with our initial sample, are illustrated in the left two top 
panels of of Figure [2j 

Based on our simulation of this survey, we find that the 
minimum heliocentric distance the original grid is sensi- 
tive to is a strong function of inclination, with d m i n ~ 22 
AU for i = 0° orbits, climbing to d m i n ~ 30 AU for 
i = 70°. This grid is in fact sensitive to inclinations as 
high as 180° outside of 36 AU. The outer edge is simi- 
larly modified by inclination, but in all cases greater than 
the 100 AU goal: at the lowest, it is sensitive to objects 
at distances as high as 164 AU for i ~ 0°, and at its 
highest it is sensitive to objects as distant as 184 AU for 
i ~ 180°. These limits are illustrated in the top right 
panels of Figure [2] 

3.2. Fraser & Kavelaars 2009 

Fraser & Kavelaars 2009 searched ~ 0.25 square de- 
grees (taken over a ~ 4 hour baseline) for TNOs with a 
grid of 95 rates and angles, aiming to be sensitive to ob- 
jects on circular orbits with heliocentric distances from 
~ 25 - 200 AU. As in Fraser et al. 2008, the adopted 



search-grid spacing was fixed in d6 and d<f>, and the final 
stacks were searched by eye. The authors state that they 
choose the (d6, dtp) spacing to limit the maximum track- 
ing error to ~ 2r; we verify that at maximum after 
4 hours, the maximum tracking error of this search-grid 
is approximately 1.25" ~ 1.8r, resulting in a maximum 
S/N degradation factor of F ~ 0.61. As the authors re- 
ported no sensitivity loss as a function of rate of motion, 
we adopt this as our uniform maximum tracking error. 

Based on our simulation of this survey, we find the 
original search-grid is over-ambitious as the high-0 re- 
gions of the grid are not populated by any real objects. 
No bound orbits observed on the ecliptic at opposition 
with heliocentric distance outside of d > 28 AU can have 
angles of motion as high as 15° from the ecliptic (Eqn. 
[ij). This, coupled with the oversampling at low-9 rela- 
tive to the maximum tracking error of 1.25" resulted in 
excess shift-vectors compared to our optimum solution. 
The optimum grid requires 28 shift- vectors, which rep- 
resents an improvement of over a factor of 3 compared 
the 95 shift-vectors searched by the authors. The orig- 
inal and optimum grids, along with our initial sample, 
are illustrated in the left two middle panels of of Figure 

m 

Similar to Fraser et al. 2008, the minimum heliocen- 
tric distance the original grid is sensitive to is a function 
of inclination, with d m i n ~ 24.5 AU for i = 0° orbits, 
climbing to d m j„ ~ 31 AU for i = 70°. This grid was 
also sensitive to inclinations as high as 180° outside of 
36.5 AU. The outer edge is significantly more distant 
than claimed: at the lowest, it is sensitive to objects at 
distances as high as 350 AU for i ~ 0° , and at its highest 
it is sensitive to objects as distant as 385 AU for i ~ 180°. 
These limits are illustrated in the middle right panels of 
Figure [2j 

3.3. Fuentes et al. 2009 

Fuentes et al. 2009 searched 0.25 square degrees (taken 
over a ~ 8 hour baseline) for TNOs with a grid of 732 
rates parallel and perpendicular to the ecliptic, aiming to 
be sensitive to objects with heliocentric distances from 
~ 20 — 200/1C/. The large number of resulting stacks 
necessitated the use of an automated detection pipeline. 
Their grid spacing was d0„ = d0± = 0".l hr^ 1 and mo- 
tions were limited to ±15° from the ecliptic. Over an ~ 8 
hour observational baseline, this grid spacing translates 
into a maximum tracking error of e ~ 0.6" ^ 0.8r, result- 
ing in a maximum S/N degradation factor of F ~ 0.76. 

Like Fraser & Kavelaars 2009, this survey also over- 
searches high-</> motions. Because of the small tracking 
error and non-optimum grid geometry, any small over- 
searched regions translate into a significant number of 
additional search vectors. The optimum grid requires 494 
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Fig. 2. — Results from survey simulations. Top panels: Fraser et al. 2008. Middle panels: Fraser &; Kavelaars 2009. Bottom panels: 
Fuentes et al. 2009. Gray distributions represent density in (da, d8) of synthetic orbits generated as described in text. Left panels include 
overlay of each survey's search grid. Middle panels include overlay of our derived "optimum" search grid. Right panels show inner and 
outer limits in heliocentric distance d vs. inclination i derived for each survey. Cross-hatched region represents limit variation due to 
eccentricity. Black triangle represents the distance at discovery and inclination of 2008 KV42- 



Table 2: Survey Characterization 



Survey 


N op t 


% Improved 




(AU) 


dmax 


(AU) 


Fraser et al. 2008 


19 


24% 


22 - 


- 36 


164 - 


- 184 


Fraser & Kavelaars 2009 


28 


71% 


24 - 


- 37 


360- 


- 390 


Fuentes et al. 2009 


494 


33% 


21 - 


- 33 


220 - 


- 245 



shift- vectors, which represents an improvement of ~33% 
over the original search grid. The original and optimum 
grids, along with our initial sample, are illustrated in the 
left two middle panels of of Figure [2j 

Similar to Fraser et al. 2008, the minimum heliocentric 
distance the original grid is sensitive to is a function of 
inclination, with d m i n ~ 21.5 AU for i = 0° orbits, climb- 
ing to dmin — 26.5 AU for i = 70°. This grid was also 
sensitive to inclinations as high as 180° outside of 33.5 
AU. The outer edge is significantly more distant than 
claimed: at the lowest, it is sensitive to objects at dis- 



tances as high as 220 AU for i ~ 0° , and at its highest it 
is sensitive to objects as distant as 245 AU for i ~ 180°. 
These limits are illustrated in the bottom right panels of 
Figure [2j 

4. EXAMPLE APPLICATION: LSST DEEP FIELDS 

A component of the LSST survey strategy will be to 
point to a single field and take hundreds of ~ 30 second 
exposures, then later combine them in software to search 
for faint moving objects (Chesley et al. 2009). On a 
single winter night, a single 9.6 square degree field at 
opposition can be observed for roughly 8 hours, resulting 
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Twotinos, 31-hr Baseline, 



-120 
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Fig. 3. — Simulated motions of near-ecliptic Twotinos after two nights (gray distribution), optimum grid for LSST observations (points) 
and non-linear corrections to grid (lines; acceleration is nearly parallel to the celestial equator for the simulated fields, but a small angle 
has been added to separate the plotted lines from each other). Left panel: Opposition observations. Right panel: Observations at 45° from 
opposition. 



Table 3: LSST Simulation Orbits 



Twotinos 




Parameter Range 


Distribution 


a 


42.8 AU 


Single- valued 


1 


25 a - 42.8 AU 


Uniform 


e 


- 0.416 


Uniform 


d 


25 a - 60.6 AU 


p(d) <x d~ 2 


i 


0° - 45° 


Uniform 



in 850 images with a combined depth of r' ~ 28. We 
will determine the feasibility of stacking multiple nights 
of data to search for specific TNO populations at even 
fainter magnitudes. 

We have run simulations of similar LSST observations 
to demonstrate the utility of our method for generat- 
ing shift- vectors. One of this method's chief advantages 
is the ability to strictly limit the orbital parameters of 
interest, being certain to sensitize the resulting stacks 
only to those orbits of interest without creating stacks at 
extraneous rates of motion. As such, the simulation de- 
scribed here is designed to detect a population that has 
some orbital parameters that are well-defined; namely, 
objects in the Neptune 2:1 resonance, or "Twotinos." 
The orbital ranges used to simulate this population are 
listed in Table 3. 

We compare observations made at opposition to those 
made at 45° away from opposition, which have increasing 
contributions from non-linear components of motion. If 
an opposition field can be observed for 8 hours on a sin- 
gle night, fields at this elongation can be observed for ~ 7 
hours per night. We simulate observing both fields for 
one, two, and three nights, resulting in the total observa- 
tional baselines listed in Table 4. The seeing is assumed 
to be the projected 75'^-percentile in r', roughly 0".89 
(Tyson et al. 2009). We adopt e max — 0.8F, similar 
to what has been used in previous automated-detection- 
based pencil-beam surveys, resulting in e max ~ 0".7. 

Table 4 contains the results of our simulations. We 



have estimated the required number of shift-vectors, 
pixel additions, and pixels to search (without the ad- 
ditional optimization from the tree of grids discussed in 
Section 2.3.4). Figure ^ illustrates the optimum non- 
linear grids generated for both elongations. Also illus- 
trated is the magnitude of the error between the final 
positions of real sources and the linear extrapolation of 
their position from their initial motion (0q,^o)- 

It is useful to compare the total required number of 
pixel additions (N a( id) and pixels to search (N aearc h) 
to the most computationally-intensive digital-tracking 
survey to date. The HST survey performed by Bern- 
stein et al. (2004) required N a dd ~ 10 16 additions and 
N searc h ~ 7 x 10 13 pixels with 2004 computer technol- 
ogy. Because the LSST PSF is significantly larger than 
that of HST, fewer shift- vectors are required for the same 
observational baseline. In our LSST Twotino simula- 
tions, even the three-night opposition solution requires 
as many or fewer N add (6.4 x 10 16 additions) and N search 
(2.5 x 10 13 pixels). Including the potential improvement 
of a factor of ~ 6.25 in the total number of pixel addi- 
tions required with the ad diti on of a second-level grid to 
the three-night arc (Eqn. |l~2"| ), it is clear that even with 
technology that will be over a decade out of date by the 
time LSST comes online, multi-night stacking of LSST 
data will be feasible for some TNO populations. 

Generalizing to the entire TNO population does not 
vastly increase the computational demand. We repeated 
the 8-hour opposition simulation allowing a to vary uni- 
formly over the same range as d and the number of 
required shift-vectors increased by approximately 15%. 
This increase is largely driven by the addition of low- 
eccentricity objects at low heliocentric distances. 

5. CONCLUSIONS 

We have shown that by generating search-grids for 
pencil-beam surveys in a probabilistic way and imple- 
menting an optimized grid geometry, significant reduc- 
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Table 4: LSST Simulation Results 



Twotinos 










/8 t b (hr) 


Depth (r' mag) 


Nshift 


Nadd 


■^V search 


0° 8 


28 


200 


5.4 x 10 14 


6.4 x 10 11 


32 


28.4 


3,359 


1.8 x 10 16 


1.1 x 10 13 


56 


28.6 


7,743 


6.4 x 10 16 


2.5 x 10 13 


45° 7 


27.9 


113 


2.7 x 10 14 


3.6 x 10 11 


31 


28.3 


1,965 


9.4 x 10 15 


6.3 x 10 12 


54 


28.5 


6,864 


4.9 x 10 16 


2.2 x 10 13 



tions (25 — 75%) can be made in the total number of 
operations required, and using a tree-of-grids structure 
can reduce the total pixel additions required even fur- 
ther. Besides improvements in computational require- 
ments, this method also ensures uniform sensitivity to 
the entire targeted volume of orbital parameter space. 

This method also provides a simple way to include 
non-linear components of motion necessary to track ob- 
jects over long observational baselines. This combination 
of reduced computational requirements and simplicity 
of extending the observational arc makes the prospect 
of stacking multiple nights of data appear completely 
feasible. The additional depth and improved precision 
of measured orbital properties derived from multi-night 
stacks would be hugely beneficial if used in upcoming 
large-area surveys like LSST and Pan-STARRS. Simula- 
tions of the LSST deep observations show that, for cer- 
tain TNO populations, it would require less computa- 
tional effort to obtain the same depth as Bernstein et 
al. (2004) over 500 times the area by stacking multiple 
nights of LSST data with the methods described here. 

Additionally, this method allows for simple character- 
ization of the orbital sensitivity of existing pencil-beam 
surveys in literature. From our characterization of the 
surveys by Fraser et al. (2008), Fuentes et al. (2009), 
and Fraser & Kavelaars (2009), we have shown that the 
orbital sensitivity of these surveys varied from what was 
advertised. Their extra sensitivity to the motions of dis- 



tant sources makes these surveys sensitive to populations 
like distant members of the Scattered Disk and Sedna- 
like objects, and their non-detection provide useful upper 
limits on these populations. We will explore these upper 
limits in future work. 

Since the discovery of extremely high-inclination, low- 
pericenter objects like 2008 KV42 (i — 110°), it is im- 
portant to understand these surveys' sensitivity to such 
populations. We note that Fuentes et al. (2009) is the 
only survey characterized here that was sensitive to ob- 
jects with the same inclination and distance at discovery 
as 2008 KV42 (see Figure [2j bottom-right panel), but we 
contend that it is unlikely that they would have recog- 
nized any detection as belonging to this retrograde pop- 
ulation. Since follow-up has been performed only rarely 
for these deep surveys, orbital inclination and heliocen- 
tric distance usually remain degenerate, with the pro- 
grade solution lying at lower distances (though this de- 
generacy is rarely acknowledged). It is possible that de- 
tections labeled as high-inclination prograde objects may 
in fact be objects with retrograde orbits like 2008 KV42 
at greater distance. Only multi-night arcs or follow-up 
at later epochs can break this degeneracy and clearly 
identify objects belonging to rare TNO sub-classes. 
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